#include <stdio.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>

#define M 2
#define N 2
#define K 3

int main() {
    cublasHandle_t handle;
    cublasCreate(&handle);
    // A: 2×3 (行主序)
    float A[6] = {1, 2, 3, 4, 5, 6}; 
    // B: 3×2 (行主序)
    float B[6] = {7, 8, 9, 10, 11, 12}; 
    // C: 2×2 结果矩阵
    float C[4] = {0};
    float alpha = 1.0f, beta = 0.0f;
    // 分配 GPU 内存
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, sizeof(A));
    cudaMalloc(&d_B, sizeof(B));
    cudaMalloc(&d_C, sizeof(C));
    cudaMemcpy(d_A, A, sizeof(A), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, sizeof(B), cudaMemcpyHostToDevice);
    cudaMemset(d_C, 0, sizeof(C));

    // A 和 B 都是行主序，因此使用 CUBLAS_OP_T
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,N,M,K,&alpha,d_B,N,d_A,K,&beta,d_C,N); 
    // 拷贝结果回主机
    cudaMemcpy(C, d_C, sizeof(C), cudaMemcpyDeviceToHost);
    // 输出结果
    printf("Result matrix C:\n");
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            printf("%f ", C[i * M + j]);
        }
    }
    printf("\n");
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    cublasDestroy(handle);

    return 0;
}

/*

cublasStatus_t cublasSgemm(
    cublasHandle_t handle,  // cuBLAS 句柄
    cublasOperation_t transa, // A 是否需要转置
    cublasOperation_t transb, // B 是否需要转置
    int m,  // 结果矩阵 C 的行数
    int n,  // 结果矩阵 C 的列数
    int k,  // A 的列数 / B 的行数
    const float *alpha,  // 标量 α
    const float *A,  // 矩阵 A 的设备指针
    int lda,  // A 的 leading dimension（行主序存储时是 A 的列数）
    const float *B,  // 矩阵 B 的设备指针
    int ldb,  // B 的 leading dimension（行主序存储时是 B 的列数）
    const float *beta,  // 标量 β
    float *C,  // 结果矩阵 C 的设备指针
    int ldc   // C 的 leading dimension（行主序存储时是 C 的列数）
);

不使用cublasSgemm transa与transb参数
由于C/C++程序中输入的A和B是按行存储，所以在的情况下，cublas其实读取到的是A和B的转置矩阵AT和BT

根据线性代数的规则可知CT = (A x B)T = BT x AT 所以cublasSgemm API中几个参数设置如下

设置了cublasSgemm的transa与transb参数=CUBLAS_OP_N
乘法左矩阵为BT=参数设置为B，乘法右矩阵为AT=参数设置为A
结果矩阵的行数为CT的行数=参数设置为N
结果矩阵的列数为CT的列数=参数设置为M
乘法左矩阵列与乘法右矩阵的行=参数设置为K
按列优先读取乘法左矩阵B的主维度（即BT有几行）=参数设置为N
按列优先读取乘法右矩阵A的主维度（即AT有几行）=参数设置为K
结果矩阵存储在参数C中，它的主维度（即有几行）= 参数设置为N
*/